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The connection between membrane inhomogeneity and the structural basis of lipid rafts has 
sparked interest in the lateral organization of model lipid bilayers of two and three components. 
In an effort to investigate anisotropic lipid distribution in mixed bilayers, a self-consistent mean- 
field theoretical model is applied to palmitoyloleoylphosphatidylcholine (POPC) - palmitoyl sphin- 
gomyelin (PSM) - Cholesterol mixtures. The compositional dependence of lateral organization in 
these mixtures is mapped onto a ternary plot. The model utilizes molecular dynamics simulations 
to estimate interaction parameters and to construct chain conformation libraries. We find that at 
some concentration ratios the bilayers separate spatially into regions of higher and lower chain order 
coinciding with areas enriched with PSM and POPC respectively. To examine the effect of the 
asymmetric chain structure of POPC on bilayer lateral inhomogeneity, we consider POPC-POPC 
interactions with and without angular dependence. Results are compared with experimental data 
and with results from a similar model for mixtures of dioleoylphosphatidylcholine (DOPC) , steroyl 
sphingomyelin, and Cholesterol. 

PACS numbers: 87.16.-b, 87.16.D-, 64.70.Ja, 81.30.Dz, 64.75+g 



I. INTRODUCTION 

The membranes of cells are composed of a variety 
of amphipathic lipids that spontaneously form bilay- 
ers in the presence of water. In nature, the lipids of 
typical plasma membranes vary by head-group, chain 
length, and chain saturation. The exact composition, as 
well as external factors such as temperature and pres- 
sure, strongly affect the properties of the membrane, 
so many experiments have been aimed at characteriz- 
ing model lipid bilayers of different compositions as a 
function of these external factors [TH28]. In particular, 
the lateral organization of ternary-component mixtures 
has been of recent interest because they contain essen- 
tial ingredients for the formation of lipid rafts, under the 
right conditions. Lipid rafts have been the subject of 
a number of reviews [29H35] where they are described 
as nanometer-scale sphingolipid and cholesterol-enriched 
domains whose importance ranges from signal transduc- 
tion to organization of bioactivity in cell membranes. 

Bilayers composed of a single lipid type can undergo a 
main-chain phase transition between a liquid disordered 
phase, l a , and a gel phase, sq, as a phase transition tem- 
perature, T m , unique to the lipid type, is traversed. The 
l a phase is characterized by highly mobile lipids with dis- 
ordered chains, while in the so state the lipids are much 
less mobile and their chains are more ordered. 

A side-by-side comparison of two bilayers of nearly 
identical composition, with only a minor difference such 
as the presence of a double bond in one of the lipid 
chains, can underscore the complex nature of real cell 
membranes. For example, dioleoylphosphatidylcholine 
(DOPC) has two mono-unsaturated 18-carbon fatty 
acid chains, while in palmitoyloleoylphosphatidylcholine 



(POPC) one of those is replaced by a 16-carbon satu- 
rated chain. As a result, the DOPC main-chain phase 
transition temperature is lower than that of the POPC 
temperature by 26° C. Additionally, the area per lipid 
of POPC bilayers, 65 A 2 [36], is significantly lower than 
that of DOPC bilayers, 72 A 2 [37 . The presence or ab- 
sence of the double bond clearly plays a role in how well 
neighboring lipids pack with each other, which in turn 
determines the physical properties of the bilayer. 

Binary mixtures of low T m phosphatidylcholines and 
high T m phosphatidylcholines or sphingomyelin (SM) can 
exhibit phase separation between the l a phase and so 
phase at temperatures intermediate to their respective 
melting points [TJ [27j [38j [39] . The phase coexistence de- 
pends upon the relative concentrations of each lipid type, 
and can be observed directly in fluorescence microscopy 
experiments (TJ [38] . 

Cholesterol (CHOL), another major membrane lipid, 
has a strong and complex effect on bilayer properties. 
CHOL has a tendency to induce order in the chains of 
Z a -phase lipid bilayers and increase mobility of so-phase 
lipid bilayers. When the CHOL concentration is above 
about 15% the resulting structural state is sometimes 
called the liquid ordered phase, Zq> an d is characterized 
by highly ordered, yet mobile chains. The phase behavior 
of binary mixtures of CHOL and other membrane lipids 
has been studied extensively [27l I4QU44] . Pan et al [45] 
find that CHOL-lipid interactions depend on the number 
of saturated chains present, affecting physical properties 
such as bending modulus, area per lipid, and order pa- 
rameter. Other studies EH EH [46] indicate that CHOL 
has a higher affinity for one lipid type over another, for 
example dipalmitoyl phosphatidylcholine (DPPC) or di- 
saturated sphingomyelin (SM) over POPC. CHOL shows 
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a unique level of specificity in inducing order because 
the two "faces" of CHOL are structurally distinct. Here, 
"face" refers to the sides of the flat fused ring structure 
of CHOL, one of which has protruding methyl groups 
(rough f3 face) , while the other does not (smooth a face) . 
Pandit et al. [47 utilize atomic-level simulations of lipid- 
CHOL mixtures to reveal that the effective molecular 
area of CHOL is smaller in POPC than in DOPC and 
even in DPPC. They attribute this effect to a combina- 
tion of the anisotropic chain structure of POPC and the 
anisotropic faces of CHOL. 

This differential behavior of CHOL in lipids must have 
a strong effect on the lateral organization of ternary com- 
ponent mixtures consisting of a high T m lipid, a low T m 
lipid, and CHOL. Experiments [H EH Ell E3 ETJ SHJ EH] 
consistently indicate that these ternary component mix- 
tures have a complex lateral structure that is highly de- 
pendent on the relative abundance of components and 
external parameters such as temperature and pressure. 
In particular, at certain concentration ratios, the liquid 
phases l a and lo are known to coexist and are the basis 
for the concept of "lipid rafts" which are believed to exist 
in biological plasma membranes [32 , 34, 35 , 50 , and are 
implicated in a number of biologically important cellular 
functions (See e.g. [3Q| 133]). 

A number of experimental treatments of ternary mix- 
tures have been published [2-23, 51 , but the details on 
the formation of rafts in POPC-SM-CHOL mixtures re- 
mains a subject of debate. A number of experimen- 
tal studies have treated specifically this system [21] [241 
[28l [50] . Veatch and Keller [1 report liquid-liquid phase 
coexistence in giant unilamellar vesicles using fluores- 
cence microscopy. Zhao et al. [15 do not observe liquid- 
liquid phase separation in POPC-SM-Chol mixtures at 
the micron level. They do observe Id and so phase co- 
existence and suggest that domains may exist in smaller 
sub-micron clusters. Overall, current experimental stud- 
ies underscore a need to carefully probe subtle atomic 
level effects in the POPC-SM-CHOL ternary system. 

Several recent theoretical models have been published 
that describe bilayer phase behavior in three-component 
systems. Putzel et al. [52] [53] have proposed a phe- 
nomenological model intended to elucidate the mecha- 
nism for phase separation in ternary lipid systems by ex- 
amining phenomenological free energy functions for mix- 
tures consisting of saturated-chain lipids, unsaturated- 
chain lipids and CHOL. They were able to reproduce 
ternary phase diagrams that are consistent with those 
proposed by experiment. Elliot et al. [54] [55] have pro- 
posed a self-consistent mean-field theoretical (SCMFT) 
model of a three-component system that treats chain in- 
teractions at an atomic level. The model employs equilib- 
rium statistical mechanical analysis to construct ternary 
diagrams with phase boundaries whose presence and lo- 
cation are modulated by interaction parameters. We re- 
cently published a model for ternary lipid mixtures that 
we applied specifically to DOPC-Steroyl sphingomyelin 
(SSM)-CHOL bilayers [56 . The model projected three- 



dimensional lipid bilayer leaflets onto two-dimensional 
fields of chain order over which CHOL could diffuse. This 
model, based on a combination of equilibrium statisti- 
cal mechanics and Langevin plus Cahn-Hilliard dynam- 
ics, described temporal organization of the mixed bilayer. 
Results were displayed on a triangle diagram that agreed 
well with experiment for DOPC-SM-CHOL mixtures. 

In this paper we present a model of ternary-component 
lipid bilayers that builds on our previous modeling 
work [43] [44] [56] [57] to include lipids with non-identical 
chains. This model differs from other computational 
and theoretical models [52-55, 58-65 in that it utilizes 
atomistic molecular dynamics (MD)-generated data as 
input to a self-consistent mean-field theoretical (SCMFT) 
model that is used to characterize the structure and tem- 
poral evolution of ternary phase diagrams. We apply 
the model to POPC - Palmitoyl sphingomyelin (PSM) 
- CHOL mixtures, but the methodology can be gener- 
alized to include any two-chain lipid ternary mixture. 
To account for the asymmetric chain structure of the 
POPC molecule, an orient ational component for POPC 
is added to the interaction energy function. We show be- 
low that the effect of including orientational dependence 
is to slightly amplify the separation of lipids into regions 
rich in PSM separated from regions rich in POPC. Com- 
parison with our SCMFT simulations of DOPC-SSM- 
CHOL [56] reveals that the degree of lateral organization 
in POPC-PSM-CHOL mixtures is reduced, and differ- 
ences between order parameters in separated domains is 
more subtle. In the next section, the POPC-PSM-CHOL 
ternary system is described. Subsequent sections present 
our results and discussion. 



II. THEORETICAL MODEL 

In our SCMFT model, g3 HQ Mi [57], a three- 
dimensional lipid bilayer leaflet is cast as a two- 
dimensional field of weighted, chain-averaged order pa- 
rameters: 

n s ~ 1 

Tig Z Z 

m=l 

where the weighting fraction ^ represents the fraction 
of dihedrals in a trans configuration (n tr ) along a single 
chain to the number of dihedrals along that chain (n s ). 
j3m is the angle between the C-H bond vector and the 
bilayer normal for carbon m for the chain at the posi- 
tion ~7\ CHOL molecules are treated as two-dimensional 
'rods' that are free to diffuse through the order parameter 
field. Coupled to this order parameter field is a composi- 
tion field whose purpose is to identify the concentration 
of lipid type at any point. Both fields are discretized 
onto an underlying lattice for computation. The overall 
model methodologies are as follows: Langevin dynamics 
are used to propagate CHOL molecules over the two di- 
mensional fields, and to locally reorient POPC molecules. 
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Cahn-Hilliard dynamics are used to model the evolution 
of lateral compositional order in the composition field. 
Mean field statistical mechanics is employed to recal- 
culate the order parameter field after each dynamical 
timestep. These methods are described in the following 
subsections. 



A. The POPC-PSM-CHOL Ternary System 

The structures of POPC and PSM are shown in Fig 
[I] POPC has two non-identical hydrocarbon chains; 
one saturated palmitoyl (16:0) chain and one unsatu- 
rated oleoyl (18:1) chain. PSM ((16:0)SM) consists of 
a palmitoyl chain attached to a sphingosine backbone. 
PSM has a higher main-chain phase transition tempera- 
ture (40° C) [66] than POPC (-3°C). For this model, we 
focus on the temperature range intermediate to the two 
phase transition temperatures, namely 30° C. 




B. Order Parameter Field and Angle Field 

The SCMFT model g3 HQ ESI 152 we previously ap- 
plied to binary and ternary mixtures of lipids with two 
chains of identical structure projected single chain order 
parameters onto a two dimensional lattice. Lattice grid 
sites were spaced at a distance calculated from MD sim- 
ulations as the average distance between two neighbor- 
ing chains. To incorporate the asymmetric structure of 
POPC molecules, we now define a 'whole-molecule' order 
parameter for each site as an average: 

Smol ~ -^(Ssnl + $sn2) (2) 

where s sn i and s sn 2 are the order parameters for the 
two chains on each molecule, as defined in Eq[l] In this 
work, we take the lattice spacing to be the average near- 
est neighbor center of mass distance between molecules. 
Although the two chains in PSM do not differ appreciably 
in order parameter value, PSM moleular order parame- 
ters are also calculated using Eq [2] to keep the model 
internally consistent. 

In order to incorporate the asymmetric structure of 
POPC molecules, a molecular orientation variable, p, is 
introduced and is associated with each lattice point that 
contains all or part of a POPC in the local concentration 
field. For a POPC molecule, p is calculated as the angle 
made between a vector which points from the atom la- 
beled 'cl' towards the atom labeled 'c2' in F ig JT] and a 
fixed direction in the plane of the bilayer (Fig[2| 






POPC PSM 



FIG. 1: The structures of POPC and PSM. The labelled ar- 
rows indicate specific atoms used to define POPC molecular 
orientation. 



In contrast to our previous implementation of the 
SCMFT model [56] for DOPC-SSM-CHOL mixtures, we 
now consider a case where the low melting point lipid, 
POPC, has nonequivalent acyl chains. To address the 
nonequivalence of POPC chains, we consider the effect of 
adding an orientation degree of freedom to each POPC 
molecule, as we describe below. Modifications to the pre- 
vious model [56] include redefining the order parameter, 
making adjustments to the Hamiltonian and mean-field 
free energy, and modifying the method of updating time 
steps. The modifications are summarized in the following 
subsections. 



FIG. 2: The angle p for a POPC molecule at site i. The 
black arrow represents the orientation vector of the lipid at 
site i, pointing from the saturated chain of POPC towards 
the unsaturated chain. We measure the angle p with respect 
to the y direction of the x — y axis of the lattice. 

The angle field is assumed to evolve dynamically ac- 
cording to a Langevin equation: 



dpi 
dt 



OF 

dpi 



(3) 



M p is the mobility of the angle as it diffuses azimthally 
at the site i and is related to the diffusion constant by 
M p = Dp/kbT. F is the free energy, shown below, and ^ 
is a stochastic noise component that satisfies fluctuation 
and dissipation relations. 
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The diffusion constant D p is found by performing sta- 
tistical analysis on measurements of the angle p cal- 
culated directly from molecular dynamics simulations. 
This angle is recalculated after a nanosecond of simu- 
lation has passed. Fig [3] shows a distribution plot of 
trajectory- averaged changes in the orientation angle after 
one nanosecond intervals for all of the POPC molecules in 
a MD simulation (details of the simulation are discussed 
below), 




FIG. 3: Distribution of the change in angle after 1 ns for 
MD-simulated POPC molecules. 

The gaussian-like shape of the change in angle over 
time allows us employ an Einstein relation to extract the 
orientational diffusion constant from MD simulations: 



(\ P (t + T)- P (t)\ 2 ) 



2D p t 



(4) 



over a short time r giving us a diffusion constant in 
units of radians 2 /ns. A best-fit line of a plot of 
(|p(t + r) — p(t)\ 2 ) vs. r, shown in Fig H yields twice 



The numerical value of D p and 



the diffusion constant. 

other parameters used in this work are given in table [I] 



Concentration Field 



To model two different types of lipids in our mixtures, 



we include a concentration field, ^i(t) 



i(t) which describes the compositional makeup of 
site i at time t. As in our earlier work, [56 , c psrnyi (t) 
is the concentration of PSM at site z, and c popCj i(t) is 
the concentration of POPC at site z, both of which range 
from to 1. Therefore, ipi(t) varies continuously from 
— 1 to +1 with the —1 representing pure POPC and +1 
representing pure PSM. Locally, ipi(t) may vary but the 
total number of lipids present in the system remains con- 
stant and so the sum over all sites, ^^(t), is fixed at 
each time. The concentration field evolves at each lattice 
point through the Cahn-Hilliard equation [67] : 



dt 



m2 dF 



(5) 




FIG. 4: Average square change in the POPC orientation angle 
as a function of time. The slope of the best fit line (dotted 
line) yields 2D P . 



F is the free energy and 7^ is a conserved stochastic noise 
component. T is the lipid mobility and is related to the 
lipid diffusion constant by T = Dn^jk^T . The value of 
Dupid is unchanged from earlier work, and is given in Ta- 
ble 1. Eq|5] ensures that, while local lipid concentration 
field values may vary in time in response to thermal fluc- 
tuations and tendency towards free energy minimization, 
the total lipid concentrations in the system are conserved. 



D. Hamiltonian 

The system Hamiltonian contributions arise from lipid- 
lipid chain interactions, lipid chain-CHOL interactions, 
and CHOL-CHOL interactions: 



Hi 



lip— lip 



-CHOL 



CHOL-CHOL- 



(6) 



The Hamiltonian couples the order parameter and con- 
centration fields to each other and to the overlying 
CHOLs. Due to the anisotropic nature of the POPC 
molecules, the first term in Eq [6] differs from our pre- 
vious model [56] and is explained in detail below. The 
second and third terms in Eq [6] have not been altered 
from [56 , but are also briefly summarized below. 



1. Lipid-Lipid, Interaction Term, Hu v -u v 

Lipid-lipid interactions are calculated following the 
model first proposed by Marcelja [68;, defined in terms 
of order-parameter based pairwise interactions between 
chains. The model was originally applied to DPPC- 
CHOL mixtures [43] [44] [56] , and produced heat capaci- 
ties and a phase diagram that agreed quite well with ex- 
periment [40]. In this application, we consider two model 
scenarios: (i) a model that contains no orientational de- 
pendence for POPC-POPC interactions, and (ii) a model 
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that includes explicit orientation dependence for POPC- 
POPC interactions. 

If we include no orientation-dependent interactions 
between neighboring POPC molecules, the lipid-lipid 
Hamiltonian for a pure POPC bilayer is written 



pope— pope 



(7) 



where the angled brackets indicate that the sum is taken 
over nearest neighbor pairs and the coupling constant, 
ypopc^ . g a ph enome nological parameter that is tuned in 
such a way that the calculated phase transition temper- 
ature of pure POPC matches the experimental value. To 
accomplish this, the average order parameter is calcu- 
lated in a system of pure POPC for a range of tempera- 
tures with a given value of Vq° pc . A plot of order param- 
eter against temperature reveals a curve that exhibits a 
sharp increase in order as temperature is decreased be- 
low a threshold temperature, T m which identifies T m as 
the phase transition temperature. Vq° pc is phenomeno- 
logically adjusted so that the simulated phase transition 
temperature is in agreement with the experimentally- 
calculated value. The numerical value for Vq° pc for this 
case is shown in table U 

With two identical saturated chains, PSM is not ex- 
pected to interact in an angular dependent fashion among 
nearest neighbors. Therefore, PSMiPSM interactions are 
written as: 



psm—psm 



(id) 



(8) 



ypsm - g ^ e coupling constant tuned in pure PSM mix- 
tures to obtain a main chain phase transition tempera- 
ture that is identical to experiment. 

For the case where we include an orientational de- 
gree of freedom for POPC molecules, we alter the POPC 
chain-chain interaction Hamiltonian as follows. We first 
recall that the concentration field variable ipi (t) can range 
continuously from —1 representing pure POPC to +1 
representing pure PSM. To incorporate the asymmetric 
nature of POPC-POPC interactions proposed above, we 
include an angular dependence to this term in propor- 
tion to the amount of POPC present between interacting 
nearest neighbors: 

^lipid-lipid = -Vo ^ S i S jl 1 + C (^i^j) X ( ( l ) ij^ji)} 
(hj) 



^2 V\ 8 ^i - v ^^j 

(ij) 



(9) 



Here, Vo is a function of fraction of POPC, Cp Q p C , and the 
fraction of PSM, c psm , present in the entire system: 



\r _ \rVopc , mm 

v — c popc ~r Vq ^sm- 



(10) 



As the fraction of POPC at site i is ^(1 — ^), the local 
normalized fraction of POPC at neighboring sites i and 
j is 



(11) 



We model the angular contribution between nearest 
neighbor POPC pairs z,j by the function X(<frij,(/)ji) of 
the relative POPC orientational angles faj and f° r 
molecules at sites i and j 



COS( 



a)- 



(12) 



The angles faj and <pji are defined in Figure [5] and dis- 
cussed in more detail in the next paragraph. The second 
term in Eq [9] couples the order parameter field to the 
concentration field and the third term couples concen- 
trations at neighboring sites. As described in [56], V\ 
and Vi are estimated directly from MD simulations. The 
concentration dependence in the first term of Eq [9] (also 
see Eq 1 1 ) ensures that interactions have an angular de- 



pendence only in proportion to the amount of POPC at 
each site. 



1 0) "ffl 



<t>ji 



FIG. 5: The angles (j>ij and <j>ji discussed in the text. Black 
arrows represent the orientations of lipids at sites i and j. <fiij 
and <j>ji are measured relative to a line connecting the lipid at 
site i with its nearest neighbor at site j. 

In order to model the condition that a lipid interacts 
more favorably with a nearest neighbor if its saturated 
chain points towards that neighbor and less favorably if 
its unsaturated chain points towards that neighbor, the 
orientation function in Eq[l2]is used. For this function, 
the orientations of POPC lipids at site i and j are il- 
lustrated with black arrows in Fig [5] which point from 
the top of the unsaturated chain towards the top of the 
saturated chain. <pij and <pji represent the orientation of 
those vectors relative to a line between site i to site j 
(the red line in Fig p| . The angle <pij is related to the 
orientation angle pi, defined above and illustrated in Fig 
|2j by a constant that depends on the location of site j 
relative to site i. If the two black arrows in Fig [5] point 
towards each other we interpret this as a saturated chain 
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of a molecule pointing towards the saturated chain of the 
neighboring molecule and faj = and =ir. This con- 
figuration is optimal and no energy penalty is imposed. 
However, if the two black arrows are pointed away from 
each other = 7r, <pj = 0), we interpret this as two un- 
saturated chains pointing towards each other, and we im- 
pose the maximum energy penalty. The constant a$ rep- 
resents the strength of the angular dependence. In this 
model, if is greater than 1, the total energy interaction 
can include a repulsion between two, nearest neighbor 
chains. Since steric repulsions are implicitly present in 
this model through mean field statistical mechanical cal- 
culations and conformation sampling, it is not necessary 
or correct to include them in the Hamiltonian. 



where sums are taken over all neighboring CHOL 
molecules [43j|44j[56]. Here, is the distance between 
two close CHOLs, with indices i and j. Ki and hvj are 
angles made between a fixed direction on the lattice and 
a fixed direction on the body coordinates of the i th and 
jth CHOLs and Kij = Ki — Kj is the difference between 
the two. From MD simulations and the success of previ- 
ous models [43[ [44] [56] , it can be surmised that a simple 
repulsive interaction is sufficient to model CHOL-CHOL 
interactions: 

t rr t tk f e sin {f^ij) ?ij ^ L 

cccc_ \0 r i:j >L 

where e = 13k B T and L is the rod length |43ll44j . 



2. Lipid-CHOL Interaction Term, Hu p -chol 



E. Mean Field Analysis 



As discussed above and described in previous publica- 
tions [43j HU [56] , lipid-CHOL interactions have an an- 
gular dependence in the plane of the bilayer because the 
"smooth" face has a tendency to induce more order in 
chains than the "rough" face does. The second term in 
the Hamiltonian accounts for this asymmetric interaction 
dependence: 

niip Tic H O L 

Hu p -chol = -J2 E VUl-AsmOi^Si (13) 



The sums are over all lipid chains and nearest neighbor 
CHOL molecules. Each CHOL 'rod' has a body coordi- 
nate system with a ^/'-axis along the length of the rod 
starting at the center of mass, and an a;"- axis extending 
from the center of mass in the direction perpendicular to 
the rod on the 'rough' side. Oij is defined as the angle 
between the ^/'-axis of CHOL i and vector connecting 
the center of mass of CHOL i and the position of lipid 
j on the lattice [44]. The coupling constant, V/ c , is found 
from from MD trajectories by linear regression analysis of 
CHOL-lipid chain interactions as a function of order pa- 
rameter, as described in previous papers [43 , 44 . The pa- 
rameter A was introduced in earlier work [44] as a means 
to incorporate the fact that lipid chains are more ener- 
getically attracted to the the smooth side of CHOL [44] . 
Thus, it serves to represent the asymmetry of lipid-CHOL 
interactions in the x-y plane of the bilayer. Numerical 
values for V\ c and A are the same as used earlier in [44] . 



3. CHOL-CHOL Interaction Term , Hchol-chol 

As in previous work, CHOL molecules are cast as two- 
dimensional rods that diffuse through the order param- 
eter and concentration fields according to a Langevin 
equation [43] [44] [56] . As in earlier work we take 

ncHOL vichol 

Hchol-chol = J2 E Kc(^)Kc(^) (14) 
i=i j=i 



After each time step, new configurations of CHOL po- 
sitions and orientations and a new concentration field 
are generated. Assuming that lipid chain order relaxes 
rapidly between timesteps, the new order parameter field 
is found in the mean-field approximation. The underly- 
ing statistical mechanical partition function is: 



E E 

i=l allconfs 



exp 



hT 



(15) 



where H tot is given in Eq[6]and described in detail above. 
Sums are over all possible configurations, represented by 
the order parameter in Eq[l] over all lipids in the system. 
It is not possible or practical to specify all single chain 
configurations, so we make use of a representative library 
of configurations relevant to a bilayer environment that 
we obtain from MD simulations [43, 44]. 

For a given set of concentration field values, {V^j^Li? 
the mean molecular field at site i due to neighboring 
lipids, CHOL molecules, and the concentration field has 
the form: 

®i = 

V 

- < sj > V [l + Ctyu^Xfa^fai)] (16) 

+ c i V lc + V 1 (C)ip i 

where v is the number of nearest neighbors to i (= 4 on 
a square lattice). Once $>i is calculated at each site, it 
is used to solve the self-consistent equations for the local 
order parameter values, within the mean-field approxi- 
mation: 



(Si) = 



^2all 



conf ^ c 



Eallconf exp\P$iS c ] 

and to find the mean-field free energy: 
F = U — TS 



E 

i=i 



k b T In Zi + H cc - - J2 ■ 



(17) 



(18) 



i=i j=i 
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F. CHOL Diffusion 

The rotational and translational diffusion of CHOL 
molecules over the order parameter and concentration 
fields are modeled by Langevin equations: 



and 



dt 



dtok 

dt 



dF 
d~? k 



duj k 



(19) 



(20) 



where rjk and are stochastic noise components and 
M r and M u are CHOL mobilities. "r^ defines the x,y 
position of the center of the k th CHOL and uj^ is the 
orientation of the CHOL body x" — y" axes introduced 
above. 77^ and ("fc are thermal fluctuations modeled as 
random variables that obey fluctuation-dissipation the- 
orems. M r and M w , the CHOL mobilities, are related 
to MD-extracted diffusion constant, D, and the rota- 
tional diffusion constant, D rotl by M r = D/ksT and 
M u = D rot /kBT [43, 44 . Numerical values for constants 
are shown in Table [U 



G. Simulation 

The system is initially given a random set of order 
parameter and concentration field values at each lattice 
site and a random distribution of CHOLs at the relevant 
concentration. The system of self-consistent equations 
is solved and the mean-field free energy is calculated. 
Following this, the angle field, concentration field, and 
CHOL molecules are updated by one time step. The self- 
consistent equations are solved again and the process is 
repeated for a total time of at least 1 microsecond. 



III. RESULTS 

In the application to POPC-PSM-CHOL mixtures, we 
consider a 100 by 100 square lattice with each lattice site 
representing a whole molecule order parameter. Molec- 
ular order parameters are averaged over both chains, so 
the field is constructed to represent 10,000 lipid molecules 
to which CHOL molecules are added. POPC:PSM con- 
centrations were simulated between either 5%to 95% or 
10% to 90%, depending on CHOL concentration, in incre- 
ments of 10%. Mixtures including CHOL were simulated 
at CHOL concentrations ranging from 0-50% in incre- 
ments of 5%. All simulated concentrations have been 
run for at least 1 microsecond at a temperature of 303 
K. To observe the effects of the addition of angular de- 
pendence, simulations were ran with values of equal 
to 0.0, 0.3, 0.5, 0.7, and 1.0. Table [E] lists the numeri- 
cal values for all input parameters in the model. While 
our goal in this modeling work was to obtain as many of 



the model parameters from MD simulations as possible, 
in some cases the mapping process involved the use of 
simulation properties (e.g. molecular energies) that were 
quite noisy. Thus the model interaction parameters that 
we used are not necessarily unique, but they are rep- 
resentative of the set of model parameters that should 
be used in this type of coarse grained modeling. In ad- 
dition to interaction parameters, the Mean Field model 
we use requires that the order parameters be symmetric 
about zero, whereas the MD library order parameter val- 
ues are not symmetrically distributed due to the bilayer 
environment. As discussed in earlier work, we introduce 



order parameter offsets 



pope 



and 



psm 



which 



are the mean values of the weighted chain order parame- 
ters found from MD simulations of pure POPC and PSM 
bilayers, respectively. These offsets are subtracted from 
each value of order parameter in the library, a necessary 
step to observe a temperature dependent phase transi- 
tion in the mean field approximation, in the absence of a 
symmetric library of chain configurations [43] [44] . 

For the MD-based parameters in Table [I] and chain con- 
figuration libraries we used the following simulation pro- 
tocols: The GROMACS simulation package [69] US was 
used to simulate a 200-lipid bilayer with compositions 
37.5% POPC, 37.5%PSM and 25% CHOL, solvated in 
7211 SPC-E water molecules. The LINCS algorithm [71] 
was used for bond constraints with an integration time 
step of 2 fs. Periodic boundary conditions were used in 
all three dimensions. Long range electrostatics were cal- 
culated using the PME algorithm [72] with a real space 
cutoff of 10 A. Van der Waals interactions were cut off 
at 16 A. The NPT ensemble was employed, allowing the 
volume to fluctuate and the Parinello-Rahman pressure 
couling scheme [73] [74] was used with a constant pres- 
sure of 1 atm. Systems were simulated for 100 ns at a 
temperature of 313 K. 

The simulated results can be qualitatively visualized 
through two dimensional color density plots of concentra- 
tion and order parameter fields. For quantitative analy- 
sis of the lateral organization we calculate binned distri- 
butions of order and concentration averaged over many 
snapshots. Examples are shown in Fig [6] The concentra- 
tion field plot for the POPC:PSM:CHOL concentration 
triplet 45:40:15 after 1 microsecond with ao = 0.0 (no 
angular dependence), Fig 6(a) , exhibits a bimodal pat- 



tern of concentration with field v alue r epresented by a 
color scheme shown in the key in |6(c)[ The order pa- 
rameter field snapshot for the same concentration after 
1 microsecond is shown in Fig |6(b)[ This field exhibits 
a pattern that is similar to the pattern of the concentra- 
tion field, with higher order parameter values associated 
with higher concentrations of PSM and lower order pa- 
rameter val ues a s socia ted with higher concentrations of 
POPC. Figs |6(e)| - |6(f)| are plots of the distribution of val- 
ues of the concentration and order fields over the lattice, 
respectively. This distribution, and others we describe 
here, was made by placing the 10,000 lattice composition 
field values into bins ranging in value from from -1 to 
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Parameter 


Value 


Method of Estimation 


Comment 




143, 160 KJ/mol 


P 


calibrated for each cxq 


-rrpsm 


63 KJ/mol 


P 




pope 


.111 


MD 


see [431144] 


psm 


.283 


MD 


see [431144] 


* 1 


20 KJ/mol 


REG 


epp ^\f\ 


Vn 


100 KT/mol 


REG 


cop t\f\ 


Vln 


3.0 KJ/mol 


REG 


see [43] [44] 


A 


2.0 


REG 


anisotropy factor see [44] 


^^oounaary 


94 nm 

\J . %J j: 11111 


MD 


Vi^rHpnrp rprjnl^i vp C1~TOTi-(j1~TOF cutoff 

11(11 V I V - V 7 1 V - 1 V - 1 J V 1 1 k J 1 V V - V 1 1 V y 1 J V 1 1 V/ 1 J \ ' V 1 I V 7 1 1 


-Lsitpta 


10 _12 777 2 /<? 

J. vy / / 1/ / o 


APPROX 

ill J_ itw^v 




D r 


10- 12 m 2 /s 


APPROX 


see [431 Hi] 




10*L> r 


APPROX 


see [Hill] 


D P 


.07 radians 2 /ns 


REG 






0.0, 1.0 


P 


dimensionless 



TABLE I: Simulation Parameters. P = Phenomenological. MD = calculated directly from MD simulations. REG = linear 
regression approximation from MD simulations. APPROX=Order of magnitude approximation from MD simulations. 
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(a)snapshot of concentration field (b)snapshot of order parameter field 



-0.6 -0.2 0.2 0.6 

(c) Concentration Field scale 



0.07 0.14 021 0.28 

(d) order parameter scale 




-1 -0.5 



0.5 1 



(e) Concentration field distribution 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

(f) Order parameter distribution 



FIG. 6: Concentration field snapshot 6(a) and scale 6(c) Order param eter s napshot 6(b) and scale 6(d) (shor t bla ck lines 
represent two-dimensional CHOL rods), concentration field distribution 6(e) and order parameter distribution 6(f) for the 
POPC:PSM:CHOL concentration 45:40:15. 



+1 (concentration) or to 0.5 (order) with width 0.01, 
for a single snapshot after lfis of simulation. How ever, 
we note that, based on both the color code for Fig|6(b)| 
and the narrowness of the distribution in Fig 6(f)| the 
variation of order over the lattice at this concentration is 
reduced compared to the variation in concentration over 



the lattice. 

To study the effect of POPC-POPC interaction angu- 
lar dependence in the model, we have run simulations 
with ao values between 0.0 and 1.0. The effect of in- 
creasing ao is to very slightly increase the magnitude 
of lateral organizational patterns in the order parameter 
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(b)a = 1.0 
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(d)a = 1.0 



0.07 0J4 0.28 

(e)order parameter scale 



035 



FIG. 7: Order parameter field snapshots for POPC:SM:CHOL concentrations 50:50:0 (top) and 40 :40:20 (bottom) for values 
of ao — 0.0 and ao = 1.0. CHOL molecules removed for clarity in bottom figures. Scale is shown in 7(e) 



field at low CHOL concentration. Fig [7] illustrates the 
effect of ap for two m ixtures, POPC:PS M:CH OL 50:50: 
(Figs [7(b)1 and [7(b)] ) and 40:40:20 (Figs [7^ and [7(d)} . 
In both cases, ao = 1.0 leads to a slightly higher de- 
gree of lateral segregation in the order parameter field. 
Generally, we find that the incorporation of angular de- 
pendence leads to minor changes in the properties of the 
mixtures for CHOL concentrations less that about 20%. 
At CHOL concentrations closer to 50%, snapshots with 
ao = 0.0 and ao = 1.0 are indistinguishable. 

By examining order parameter and concentration field 
distributions, and snapshots, we can locate points on a 
ternary triangle diagram where the model predicts bi- 
modal distributions of order and concentration, which 
we then relate to experimental ternary bilayer triangle 
diagram regions. Fig[8^a) is a ternary plot published by 
Veatch and Keller [1 , based on fluorescence microscopy 
for the temperature of 23° C. Superimposed over their 
plot are arrows linking simulated concentration points 
on the diagram to system snapshots in figures |8jb)-[8^k), 
for simulations carried out with ao = 1. In comparing 
our model to experimental data based on fluorescence 
microscopy, it is important to note that the simulation 
time scale, on the order of 1 /is, is much smaller than ex- 
perimental time scales. Lateral inhomogeneities in this 



model as well as those in other models based on the Cahn- 
Hilliard [67j [75] equation, tend to increase in size over 
time under favorable conditions. 

Fig[8^b) shows the final order parameter snapshot for 
the POPC:PSM:CHOL concentration 30:50:20 and Fig 
|8jc) shows the final concentration field snapshot the same 
concentration. This concentration corresponds to the 
point labeled (b,c) with an arrow in Fig[8^a). We can see 
in the concentration field rounded regions of light green 
color, indicating they are enriched with POPC. They are 
embedded in a background of a dark blue color, indicat- 
ing enrichment in PSM. The order parameter field color 
plot for this mixture, Fig[8jb), has a pattern that is sim- 
ilar to the concentration field plot, but note that the 
difference in the magnitude of the order parameters, the 
color scale difference, is reduced. 

On the triangle diagram, experimental data points sur- 
rounding the point (b,c) are black circles, which denote 
coexisting high and low order domains. In the ternary 
plot of the same POPC-PSM-Chol system published by 
de Almeida et al. [27], this concentration point is con- 
sidered to be within a three-phase region of coexisting 
gel, liquid disordered and liquid ordered phases. Zhao et 
al. [10], argue that in some cases micron-sized separated 
regions may be light induced [10 , but they do not rule 
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(I) Order parameter scale ( m ) Concentration field scale 



FIG. 8: Experimentally derived phase plot [T, superimposed with labels marking the positions of the POPC: PSM:CHOL 
concentration fields and order parameter fields discussed in the text. In the simulation snapshots, CHOL molecules are 
removed for clarity. In all cases, cto = 1.0. 



out their existence on a smaller size scale. 

In Fig (8^d) we show the final order parameter field 
snapshot and in Fig [8je) we show the final concentra- 
tion field snapshot for the 60:30:10 mixture. The corre- 
sponding point on the triangle diagram is labeled (d,e) 
by an arrow. Again, we see the formation of higher or- 
dered regions in the order parameter field coinciding with 
PSM-rich regions in the concentration field and regions 
of lower order coinciding with POPC-rich regions. The 
experimental data points on the triangle diagram in Fig 
[8^a) shows black dots at regions surrounding this con- 
centration point, representing the coexistence of liquid 
ordered and liquid disordered phases. 

In Figs^f) and|8jg), we show the final order param- 
eter and concentration field snapshots at the concentra- 
tion triplet 60:20:20, labeled (f,g) with an arrow on the 
triangle diagram. For this point, where the POPC con- 
centration is high, there are no clear bimodal patterns in 
either field. However, CHOL does have an effect on the 
order, raising the values of chain order without forming 
large-scale lateral inhomogeneity. At similar POPC:PSM 



ratios (3:1), but in the absence of CHOL (not shown in 
this figure) we find that the order parameter field is quite 
uniform with an average value of ~0.15, with a very nar- 
row distribution. However, Figs[8jf) and|8jg) show that 
20 % CHOL increases the order into the range ~0.16- 
0.17 with a wider distribution. In the order parameter 
field color density plot, there are diffuse regions that are 
more heavily populated by higher order and diffuse re- 
gions that are more heavily populated by lower order, but 
sharp demarcation between regions is absent. The ex- 
perimental data points near point (f,g) are white circles 
that represent a single liquid phase region. De Almeida 
et al [27 and Pokorny et al. [24] find evidence for coex- 
istence of liquid ordered and liquid disordered phases at 
this point on their POPC-PSM-CHOL triangle diagram. 

The concentration triplet 10:50:40 is shown in Fig[8pi) 
and[8ji) and labeled (h,i) on the triangle diagram. At this 
concentration there are no significant inhomogeneities in 
either the order parameter or the concentration fields. 
In this case there is a greater presence of PSM and the 
result is a uniformly higher order parameter field and a 
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uniform concentration field. The presence of CHOL at 
this point has the effect of inducing order on all of the 
chains. The experimental data points surrounding the 
simulated concentration are white circles that represent 
the presence of a single liquid. Pokorny et al. [24] and 
de Almeida et al. [27] also identify this region as a single 
liquid disordered phase. 

Figs|8jj) and 8^k) show order parameter and concen- 
tration field final snapshots for 80% PSM and 20% POPC 
in the absence of CHOL. The corresponding concentra- 
tion is labeled (j,k) on the experimental diagram. There 
is no discernible separation of order parameter and con- 
centration field values between high and low values. Both 
order parameter and concentration field distributions for 
this ratio are quite uniform. At this same concentration, 
the experimental data point shown in Fig [8^ a) is a white 
square, indicative of a single highly ordered gel phase. In 
our simulations, along the 0% CHOL concentration axis, 
as POPC concentration increases from 10% to 90%, the 
order parameter field evolves smoothly (within the reso- 
lution the number of simulations) from a high value, to 
low value with increasing POPC. Bimodalities in order 
are seen as shoulders rather than separated peaks, and 
one can discern subtle variations in snapshots, as seen 
for example in Fig [7] In the concentration field, however, 
there is a region of coexistence between PSM concen- 
trations of 0.30 and 0.75 at 0% CHOL. At 50:50:0, the 
domains in the model are no longer circular, but perco- 
late across the simulation box (see Fig [7]). The triangle 
diagram published by de Almeida et aL\21\ exhibits bi- 
modal behavior along the 0% CHOL concentration line, 
with a gel phase at the high PSM end, liquid disordered 
phase at the high POPC line, and a broad coexistence 
region between. 



CHOL 



0.35 



POPC 




PSM 



FIG. 9: POPC-PSM-CHOL ternary phase diagram with a = 
1.0. The shaded region indicates the presence of bimodality in 
the order parameter distributions. Dots indicate the average 
order parameter values, with color key to the right. Two 
overlapping dots that appear in the bimodal region represent 
the average order parameter of each peak. 



In Fig [9] we collect all of the results of 1 {is simula- 
tions onto a triangle diagram based on order parameters 
for cyq = 1.0. The triangle diagrams for other values of cyq 
are not qualitatively different. The colored dots represent 
the average value of molecular order parameters, with the 
color of the dots denoting the magnitudes. Within the 
shaded region, there are two dots at each simulated mix- 
ture point representing the respective values at each peak 
in order parameter distributions which show bimodality. 
The bimodal region itself is outlined by a boundary con- 
tour that was approximated by fitting distribution peaks 
to gaussians and interpolating where single peaks would 
branch into double peaks. The resolution of this bound- 
ary is limited by the number of simulations carried out. 
As expected, Fig [9] shows that increasing CHOL concen- 
tration has the effect of increasing the order parameter 
for the full system. 

As can be seen in the ternary figure, the patterns of 
the colored dots exhibit a shift from high order parameter 
value (orange, yellow) in PSM-rich regions to low order 
parameter (green, blue) in POPC-rich regions. At con- 
centrations of CHOL below about 0.3, the order param- 
eter and concentration fields separate into regions rich in 
POPC and regions rich in PSM, as was discussed above. 
To quantitatively analyze this behavior, we consider the 
15% CHOL concentration line in the ternary plot, whose 
distributions of or der paramet er and concentration fields 
are shown in Figs |10(a)ll0(p)| At low POPC concentra- 
tion, 5:80:15 ( Figs 1 1 ( a) | and 1 1 (b )] ) , the order fields have 
a single peak and the average value occurs at approxi- 
mately 0.28 (a yellow dot in Fig|9|. As we move toward 
higher POPC concentrations, the order param eter va lue 
is redu ced. At the concentration 35:60:15 (Figs 10(c)] and 
|10(d)[ ), we start to see the emergence of a second peak 
in concentration around 0.2. Moving further into the 
bimodal region, this concentration field peak grows and 
event ually d omi nates the distribution. Above 55:30:15 
(Figs [l0(k)1 and (10(1)1 ), again, only a single peak is visi- 
ble in the distribution. Note that in all cases the order 
parameter field exhibits a much more subtle level of bi- 
modality, through shoulders in a single peak rather than 
separated peaks. The addition of angular dependence 
does not change the distributions significantly. 

As CHOL increases above 15% we find that the order 
parameter bimodal peaks begin to merge, while the con- 
centration field distribution remains bimodal. The dis- 
tributions of the order and concentration fields at 50% 
CHOL are shown in Fig [TT] At high PSM concentra- 
tions, 10:40:50 (Figs |ll(a)| and |ll(b)[ ) there is a single 
peak in both the concentration field and the order pa- 
rameter field. This is also true at high POPC concentra- 
tions, 40:10:50 (Figs 1 l(g)| and |1 1(h) > ). How ever, be tween 
the se two extremes, (Figs |1 1(c) and |ll(d) and Figs ll(e)| 
and ll(f)[ ), the concentration field distributions split into 
two clearly separate peaks, but the order parameter dis- 
tributions still show just a single peak. In this model, 
we find that the presence of CHOL at moderate to high 
concentration (above 30% CHOL concentration) has the 
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(f)25:60:15 concentration 
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(o)75:10:15 order 




(p)75:10:15 concentration 



FIG. 10: Distribution plots of order 
labeled for the last snapshot of 1/is 



parameter and concentration field for the POPC:PSM:CHOL concentrations, for ao — 1.0, 
simulations along the 15% CHOL concentration line. 



effect of ordering lipid chains of both PSM and POPC to 
such a degree that we are unable to distinguish between 
POPC and PSM by chain order. 



IV. DISCUSSION 

We have presented a quantitative model for ternary 
mixtures of POPC, PSM, and CHOL. The model rep- 
resents a high level of coarse graining that allows us to 
simulate many different mixtures on relatively large lat- 
tices for relatively long times. The results presented in 
the form of distributions of order and concentrations, 
and in the form of final snapshots, fit well with exist- 
ing experimental data. The degree of phenomenology 
is necessarily greater than our earlier work [56] but ev- 
ery effort was made to use input from MD simulations 
where possible. As we found in earlier work, the order 



parameter field and concentration field distributions vary 
according to the mean field approximation, and are ex- 
pected to provide an accurate description of the system 
as long as one is not near a critical point. In the model, 
as relative concentrations of POPC, PSM, and CHOL are 
changed, bimodalities in the concentration distributions 
can be found over a portion of the triangle diagram as 
illustrated in Fig [9] For intermediate CHOL concentra- 
tions, bimodalities in order parameter field distributions 
tend to parallel those found in concentration field dis- 
tributions. However, in the POPC-PSM-CHOL results, 
differences in chain order between adjacent domains are 
generally smaller in magnitude when compared to results 
from the same model applied to DOPC-SSM-CHOL [56]. 
At higher CHOL concentrations (above about 30%) or- 
der parameter bimodality is not found even for mixtures 
that exhibit concentration field bimodality. 

A key difference between the results of this paper and 
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(a)10:40:50 order 
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(b) 10:40:50 concentration 
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is practically zero. Therefore, a prediction of this model 
is that molecular chain order differences may be quite 
subtle and even undetectable whereas concentration field 
distributions still exhibit bimodality, or phase separation. 



We note that, while the inclusion of an orient ational 
interaction for asymmetrical POPC molecules is a reason- 
able addition to the model, the results in our work are 
only weakly dependent on this interaction. The inter- 
action function that included POPC orientation, Eq. [9j 
was chosen as a simple way to express this interaction. 
It is possible that a different choice for the orientational 
energy expression could produce different results when 
the concentration of POPC is large. 



(c)20:30:50 order 



(d) 20:30:50 concentration 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



(g)40:10:50 order 




(f)30:20:50 concentration 




-1 -0.5 0.5 



(h)40:10:50 concentration 



In summary, our extended model for ternary compo- 
nent lipid bilayers of POPC-PSM- CHOL shows good 
agreement with experiment, and provides possible new 
insights into the atomic level compositional and order 
distributions. This SCMFT model can be utilized to de- 
scribe any system consisting of two long-chain lipids in a 
ternary mixture with CHOL, independent of chain type, 
length, or level of unsaturation, while also taking into 
account any angular dependencies. 



FIG. 11: Distribution plots of order parameter and concen- 
tration field for the POPC:PSM:CHOL concentrations, for 
ao = 1.0, labeled for the last snapshot of 1/js simulations 
along the 50% CHOL concentration line. 



our previous implementation of the SCMFT model for 
DOPC-SSM-CHOL mixtures [56 is the reduced area oc- 
cupied by the bimodal regions the POPC-PSM-CHOL 
ternary plot (see [56 ), compared to the bimodal regions 
in the DOPC-SSM-CHOL ternary plots. In the DOPC 
system, the bimodal region is wider, implying that the 
DOPC system can more readily separate in to high/low 
order regions and the inclusion of a saturated chain in 
POPC reduces the size of the regions over which sep- 
arated domains can form. Fig [12] illustrates this com- 
parative difference 60:30:10 mixtures (Fig 12(a)). In the 
DOPC mixture, small rounded domains have formed that 
are rich in ordered SSM and are surrounded by a DOPC 
background with a lower order parameter value. In the 
POPC mixtures, the color density plot shows that the 
difference in the order parameters between the regions 
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(a)DOPC 60:30:10 



(b)POPC 60:30:10 a = 0.0 



(c)POPC 60:30:10 a = 1.0 
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FIG. 12: Order parameter field final snapshots for DOPC:SM:CHOL 60:30:10 and POPC:PSM:CHOL mixtures at the same 
concentration ratios with ao = 0.0 and ao = 1.0. Scale shown in 12(d) 
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